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a  scheme  for  downdating  the  Szego  polynomials  and  given  le'  *t  squares  approximant  when  a 
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Ilessenberg  matrices.  We  describe  a  data-fitting  application  that  illustrates  how  our  scheme  can 
be  combined  with  the  fast  Fourier  transform  algorithm  when  the  given  nodes  are  not  equidistant. 
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Downclating  of  Szego  Polynomials 
and  Data  Fitting  Applications 

Ci.S.  Amni.d"  \V,B.  Gragg^  L.  lleicheB 


Abstract 

Many  algorithms  for  polynomial  least  squares  approximation  of  a  reaUvalued  function 
on  a  real  interval  determine  polynomials  that  are  orthogonal  with  respect  to  a  suitable  inner 
product  defined  on  this  interval.  Analogously,  it  is  convenient  to  compute  Szego  polynomials, 
i.e.,  polynoini.nls  that  are  orthogonal  with  respect  to  an  inner  product  on  the  unit  circle, 
when  appro.xiiuaiing  a  complex<valued  function  on  the  unit  circle  in  i,he  least  squares  sense. 
It  may  also  be  appropriate  to  determine  Szego  polynomials  in  algorithms  for  least  squares 
approximation  of  real-valued  periodic  functions  by  trigonometric  polynomials.  This  paper 
is  concerned  with  Szego  polynomials  that  are  defined  by  a  discrete  inner  product  on  the 
unit  circle.  present  a  scheme  for  downdaiing  the  Szego  polynomials  and  given  least 
squares  npproxiinant  when  a  node  i.«  doletc<l  from  the  inner  product.  Our  scheme  uses  the 
Qll  algorithm  for  uniliiry  upper  llesseaberg  matrices.  We  describe  a  data-fitting  application 
titat  illustrates  how  our  sclieme  can  be  combined  with  tin-  fast  Fourier  transform  algorithm 
when  the  given  nodes  arc  not  equidistant.  Application  to  sliding  windows  is  discussed  also. 


1  Introduction 

l.et  distinct  nodes  on  the  unit  circle  and  h't  positive 

weights.  Introduce  for  comple.x-valued  functions  g  and  h  defined  .ii  ilie  nodes  the  discrete  inner 
product  on  tlio  unit  circle 

I/I 

(uJi)m  :=  '^'A:kM=k}tvl,  (1.1) 

whore  the  bar  dcnotc.s  complex  conjugation.  Polynomials  that  are  orthogonal  with  respect  to 
an  inner  product  on  the  unit  ciiclc  are  known  as  Sugo  iKthjnomials,  Let  denote  the 

family  of  orthonnrmal  Szego  polynomi.ils  with  ic.s|)oct  to  the  inner  product  (1.1),  \\hore  dj  is  of 
degree  j  and  has  a  positive  leading  coefficient.  'I'he  polynomials  <i>j  satisfy  the  Szego  incurrence 
rrinlions 


Oo(;)  =  <3o(c)  =  1/fTo, 
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orj+i<^j+i(c)  =  :Oj{z)  +  -,j+iQj(z),  0<j  <m-l,  (1.2) 

CTj+iQj^iiz)  =  :-,j+i^j(:)  +  6j{z), 

where  the  recurrence  coefficients  7^+1  €  C  iiml  cr,+i  >  0  are  determined  by 


7,+,  =  -d. (1.3) 
(Tj+x  =  (l  -  l7;+il')'^*-  ()<j<m. 

Sec,  for  c.\am]>lc,  Grenunder  and  Szoj^o  [  i.  Chapter  2],  It  can  be  shown  by  induction  that  the 
au.xiliary  polynomials  <t>j  in  (1.2)  satisfy 

dj(c)  =  ()<j<m, 

where  ij{z)  denotes  tl»e  i)olynomial  ol)taino(l  by  conjugating  the  coefficients  of  <,‘>j{z)  in  the 
power  basis.  Since  the  measure  on  tlte  unit  circle  that  defines  (1.1)  has  precisely  m  points  of 
increase,  we  have  |7;|  <  I  for  I  <  j  <  m  and  |7,„|  s  1.  The  coefficients  7^  are  hnown  as  Scliur 
immmeters,  and  we  refer  to  the  <Tj  as  the  associated  complementary  parameters,  .\lthough  the 
complementary  parameters  arc  mathematically  rodumlant,  we  retain  them  during;  computations 
in  order  to  avoid  numerical  instability.  .Vote  from  (1.3)  th.at  nl  is  the  total  weight  of  the  measure 
that  defines  (1.1),  and  that  is  the  leading  coelficient  of  for  0  <  j  <  m. 

Tor  later  reference  wo  also  define  the  polynomial 

( ■•  )  ’  ~  *  Om  - 1  (  -  )  d"  7m  4^m  — 1(^)"  (Cd) 

Then 

m 

=  n<* (I'S) 

In  particular,  -  0  fo>‘  0  <  j  <  »i- 

Example  1.1.  Let  zi,  :=  e.\p(2T;(A:  -  l)/7?i)  and  w(  :=  1  for  I  <  k  <  m,  where  i  :=  \/^. 
Then  (Tq  =  =  1  and  7^  =  0  for  I  <  j  <  m,  and  7„,  =  1.  Thus.  (i>j{z)  =  for 

()<  j  <  m,  and  =  7n”*/-(;”*  -  1).  □ 


:=  ifj. 


ft) 


1/2 

>n  ’ 


Introduce  the  discrete  iionn 


and  let  fln-i  denote  the  set  of  all  polynoiniaU  of  <le”ree  ,at  most  n  -  1.  Let  /  be  argiven 
complex- valued  function  defined  at  the  nodes  zt;.  and  consider  the  problem  of  computing  the 
polynomial  p„_i  €  fln-i.  for  some  «  <  m.  that  satisfies 

||/-/>n-ilU  =  inin  ll/-/>l|m.  (1.6) 

The  solution  Pn-i  of  the  minimization  problem  ( l.b)  can  be  expressed  conveniently  in  terms  of 
the  Szego  polynomials  Oj  as  follows.  Introduce  the  vector 

f  =  (wi/(C|).(r2/(r..) . 1  (1.7) 

and  the  m  x.  ti  iiiairix  Q  =  (</*•,, )i 

t/kj  :=  I  <  <  'o.  1  <  ji  <  «.  (1.8) 

where  u'l-  :=  y'lcj.  Note  tliat  tlie  matrix  Q  lia>  orilioiiormal  columns,  i.e..  Q‘Q  -  1.  where 
O'  ;=  ,|((j  vector  a  =  [ni.u;> . m,,)'  bo  defined  by 

a-.-Q't.  (1.9) 

Tlicii  the  polynoir.inl  p„_i  can  be  written  as 

p/l-l  ~  t .  tl.lO) 

,=i 

where  the  Fourier  coellicients  arc  independent  nf  a. 

It  is  the  purpose  of  this  paper  to  present  an  algorithm  for  downdatint)  the  recurrence  coef¬ 
ficients  7j  and  Oj.  as  well  as  the  Fourier  coellicieiiis  o,.  when  an  node-weight  pair  is  removed 
from  the  inner  product  ( 1.1).  Our  algorithm  is  ba.sed  on  the  observation  that  the  columns  of  the 
matrix  Q  arc  eigenvectors  of  a  unitary  llc.ss('iib('ig  matrix  defined  by  the  recurrence  relations 
(l.2)-(1.3).  This  makes  it  ,'o.ssiblo  to  downdale  ilie  coellicients  7;,  Oj  and  (Xj  by  applying  the 
()ll  algorithm  for  unitary  Fessenberg  matrices..  piavseiiW'd  in  [5],  with  the  node  to  be  removed 
as  shift.  Details  are  described  in  Section  2. 

\Vc  remark  that  the  problem  of  )i/«/«/mj/  ihe  (oellicient.s  7^,  Oj  and  aj  when  a  node- weight 
pair  {^m+i .  ‘t^m+i }  added  to  the  inner  product  ( 1.1 )  is  discussed  in  [8].  The  updating  problem 
can  1)0  solved  by  using  an  inverse  Qll  algorithm  for  itnitary  Hessenberg  matrices:  see  [8,  Ij. 

.\ssumc  that  wo  wish  to  (Ictcrmino  the  polynomial  /)„_i,  given  by  (1.9)-(1.10)  with  n  <  m 
when  the  set  of  m  nodes  in  the  inner  product  (l.l)  is  a  subset  of  the  set  of  eciuidistant 


nodes  {exp(2rij7.V)}j\.o‘,  with  k  :=  :V  -  m  >  0  small.  The  weigiits  wl  arc  tdl  assumed  to  be 
unity.  Then  it  may  be  attractive  to  compute  In-  (irsi  computing  the  polynomial  interpolant 
/>.v_i  €  n,v-i  on  the  sot  of  .V  o(|iiiilistimt  nodes  by  lining  the  fast  Fourier  transform  (FFT) 
f'dgorithm,  and  then  determining  from  p.v_i  by  applying  our  downdating  scheme.  The 
Fourier  coefficients  of  the  polynomial  appro.ximant  />„_i  are  then  equal  to  the  first  n  Fourier 
coefficients  of  Pm-i-  This  approach  can  similarly  be  applied  to  trigonometric  polynomials. 
Details  are  described  in  Section  3.  Computed  c.xampies  are  presented  in  Section  4  and  Section  5 
contmns  a  summary. 

Updating  and  downdating  of  polynomials  appro.Nimants  pn-i  when  all  the  nodes  tk  are 
i-eal  has  received  considerable  attention  in  the  literal nro:  see  Scott  and  Scott  [9]  and  references 
there.  A  collection  of  algorithms  for  updating  and  downdating  based  on  orthogonal  polynomials 
is  presented  by  Flliay  et  al.  [3].  .Mgorithms  for  updating  are  also  considered  in  [G,  7]. 


2  Downdating  ofSzegd  polynomials 

The  connection  between  the  Szogo  polynonuals  detennined  by  (l.l)  and  a  unitary  llesseU' 
berg  matri.x  can  be  seen  as  follows.  Using  the  basis  of  Szegd  polynomials,  we  can  write 

I  <;' <  (2.1) 

la\ 

where  for  1  <  j  <  /».  and  =  1..  hot  //*,  :=  0  for  3  <  ;+2  <  A:  <  m,  and  define 

the  upper  llcsscnborg  matrix  II  =  <l<.’liuo  ilic  unitary  matrix  U  - 

//*.,:=  0,-1  .  (2.2) 

Substitution  of  c  =  Zk,  1  <  A;  <  in.  into  (2.1)  yields  the  o(iuation 

Il'^U'^  =  (2.3) 


or,  equivalently, 


II  =  U'XU, 


(2.4) 


where  A  :=  di.igfci. c.. . c,„].  Thus.  II  is  a  uiiiiary  u|)|)or  llcssenbcrg  matrix  with  positive 

subdiagonal  elcmonts  that  is  uniquely  dotormiiK'd  by  ihe  inner  product  (1.1).  .\lso  note  that 
the  first  n  columns  of  U  make  up  the  matrix  0  th’finod  hv  ( l.S),  Algorithms  for  the  solution  of 
the  least-squares  |)iol)lem  ( l.G)  can  thei('fore  lie  \i(>\\('<l  in  l('|•ms  of  the  spectral  decomposition 
(2.-1), 


-i 


It  is  fairly  straightforward  to  sliow.  by  using  tiie  recurrence  relations  (1.2)-(1.3),  that  H 
can  be  written  as  a'product  of  in  elementary  unitary  transformations  that  are  defined  by  the 
recursion  coefficionis  and  for  the  Oj.  Wo  liavo 

=  C*  i(';i  )62(*2)  ■  ■  (2.5) 

where  the  1  <  j  <  are  m  x  m  Givens  matrices 

h-x 

-li 

fm-j-l 

and  !•'*  the  diagonal  matrix 

diag[i.l . i.-7m]. 

W'e  refer  to  the  reiuesentation  (2.5)  ns  tlie  Scliur  parametric  form  of  H. 

The  development  of  ofTicicnt  algoriilims  for  cigenproblems  for  unitary  Hessenberg  matrices  is 
facilitated  liy  the  fact  that  every  unitary  upper  Hessenberg  matri.>{  with  nonnegative  subdiagonal 
elements  has  a  uniipic  Schur  parameterization.;  I  nr  example,  when  the  implicitly  shifted  QR 
algorithm  is  applied  to  find  tlic  spectral  decomposition  of  a  unitary  Hessenberg  matrix,  a 
sequence  of  intermediate  unitary  Hessenberg  mi\trices  is  generated  that  converges  to  a  diagonal 
matrix.  In  (.5)  the  QR  algorithm  for  unitary  Hessenberg  matrices  is  formulated  in  terms  of  the 
Schur  parameters  of  the  intermediate  matrices.  This  results  in  an  implementation  that  requires 
only  0(m)  arithmetic  operations  per  iteration  on  a  matrix  of  order  m. 

Assume  for  the  moment  that  the  matrix  II  and  scaling  factor  (Tq  are  given,  but  that  the 
nodes  Zk  and  weights  lul  of  the  inner  product  (1.1)  are  not  explicitly  known.  Then  it  follows 
from  (2.4)  and  (2.2)  that  the  nodes  Zk  are  the  eigenvalues  of  II,  while  the  normalized  weight 
Wk/<To  is  equal  to  the  first  component  of  the  normalized  eigenvector  corresponding  with  z*.  (We 
.issiime  that  each  (-iiieiivector  is  scah’tl  lo  have  niiii  Kuclidean  length  and  a  nonnegalivc  first 
component.)  The  uidtary  Hessenberg  ()R  algoiiihm  can  be  used  to  compute  the  matrix  A  and 
vector  L'ei.  withoiii  computing  the  le^l  of/'.  v\iih  noi  more  than  0(m)  arithmetic  operations 
per  iieralion.  Thioi.ghont  this  p.ipei  Cj  denoi(".  ilu'  /'''  .ixi...  vector  in  C”‘.  The  QR  algorithm 
dcicrmines  one  node- weight  pair  at  a  lime,  .and  for  each  pair  computed  the  order  of  the  unitary 


0 


upper  Hessenberg  matrix  is  reduced  by  cue.  so  that  the  reduced  Hessenberg  matri::  corresponds 
with  the  node-weight  pairs  that  have  not  yet  been  determined. 

We  remark  that  in  the  case  that  the  nodes  of  the  discrete  inner  product  are  real,  then  the 
analogue  of  the  matrix  11  is  a  real  symmetric  tridiagonal  matrix  T  with  positive  subdiagonal 
elements.  This  matrix  T  contains  recurrence  coefTicients  for  orthonormal  polynomials  that 
satisfy  a  three-term  recurrence  relation.  Golub  and  Welsch  [4]  proposed  the  use  of  the  QR 
algorithm  for  symmetric  tridiagonal  matrices  for  tliu  computation  of  the  node-weight  pairs 
associated  with  T.  This  algorithm  also  determines  one  node-weight  pair  at  a  time,  and  reduces 
the  order  of  T  when  such  a  pair  has  been  found. 

Conversely,  the  construction  of  II  from  the  inner  product  (1.1)  can  be  regarded  as  an 
inverse  eigenvalue  problem.  In  particular,  given  the  matrix  A  =  diag[ji,C2,. . . ,c,nl  tmd  the 
vector  qo  :=  <^0 '[tt’i,  tvj, ....  "’e  can  perform  a  sequence  of  elementary  unitary  similarity 

transformations  whose  composition  results  in  an  in  x  m  unitary  matrix  U  such  that  the  matrix 
on  the  right-hand  of 


is  of  upper  Hessenberg  form  with  positive  suhdiagonal  elomcnts.  (The  ★  rcprose  its  an  arbitrary 
scalar  that  remains  unchanged.)  In  other  words,  V'<\q  «  rroei,  and  U'MJ  is  a  unitary  upper 
Hessenberg  matrix  II..  Consequently,  II  is  the  maiiix  corresponding  with  the  inner  product 
(l.l). 

The  transformation  of  A  to  II  can  be  pi  ifumed  using  0{rn})  arithmetic  operations  with 
the  inverse  unitary  ()ll  (lUQR)  algorithm  dob(  iil)t’<l  in  [1].,  The  lUQR  algorithm  is  an  updating 
procedure  because  it  incorporates  nodc-wciglit  pairs  one  at  a  time.  After  the  stage  of  the 
algorithm  has  boon  e.\e(  uted,  the  j  x  j  unitary  llossonborg  matrix  corresponding  with  the  inner 
product  determined  by  the  first  j  nodes  and  weights  lias  been  obtained. 

In  [8]  the  appro.simation  problem  (1.0)  is  solved  using  the  lUQR  algorithm  to  construct 
the  Schur  parameters  of  the  unitary  lles.senbcrg  matrix  11.  The  elementary  unitary  matrices 
are  accumulated  against  the  vector  f  during  the  algorithm  to  obtain  the  vector  uf  Fourier 
coelHcieuts  a  =  U'f  without  explicitly  forming  U.  in  this  way  we  obtain  the  interpolating 
polynomial  pm-i  thttt  is  the  solution  of  (1.6)  with  u  =  m  in  0{m'^)  arithmetic  operations.  Of 
course,  the  partial  sums  of  the  Fourier  expansion  oi  the  interpolating  polynomial  yields  the 
solution  (1.10)  of  (l.G)  for  each  ;t  <  in.  Moreover,  if  one  is  interested  in  computing  only  p„_i 

0 


for  n  <  m.  then  the  algorithm  can  be  curtailed  that  only  0{mn)  arithmetic  operations  are 
required  for  the  computation  of  the  parameters  See  [8]  for 

derails. 

Xow  assume  that  wo  have  solved  the  least->«|uares  problem  (1.6)  with  n  =  m  by  the  method 
described  in  [8|.  so  that  .sets  of  Sdmr  parameters  and  of  complementary  parameters 

corresponding  with  the  inner  product  (1.1).  as  well  as  sets  of  Fourier  coefficients 
of  the  interpolating:  polynomial  pn,.|  arc  c.'cplicitly  known.  In  the  downdating  problem,  we  seek 
to  solve  the  corresponding  least-s<|uares  problem  wit  h  one  term  deleted  from  ( 1.1 ).  In  particular, 
let  r  be  the  node  tliat  is  to  be  ileleted  from  the  inner  product  (1.1)  and  let  ib  be  the  sqnare*root 
of  the  associated  weight..  In  order  to  .’.implify  some  formulas  that  follow,  we  assume  without 
loss  of  generality  that  i  =  and  w  -  ir„,. 

Introduce  the  inner  product 

■-  ^ rju'f.  (2.7) 

t  si 

and  discrete  norm 

!!vllm-i  :=  (//..7)m!l- 

Wo  seek  the  polynomial  €  llm-i  such  that 

=  lldll  ||/-Hlm-J-  (2.8) 

a€tl..i_a 

Denote  the  fatnily  of  orthonormal  Szego  polynomials  with  positive  leading  cocfricient  associatod 
with  (2.7)  by  {Oj}"'jJ.  .Vlso  let  (Tj}”!,"'  and  {'t,  denote  the  sets  of  recurrence  coefficients 

for  the  0j,  and  let  ,\  ;=  diag[ci . and  <|„  :=  'T,7'[n'i . where  (t,,  :=  (ffj^  - 

is  the  total  weight  of  the  measure  that  defines  (2.7).  Then  from  the  above  discussion, 
there  is  a  unique  unitary  l(cssonl)org  matri.K  //  and  a  unique  unitary  matrix  U  such  that 

Uqo  =  <»oet 


and 


II  =  U'\U, 


(2.9) 


where  donates  the  axis  vector  in  .Moreover,  the  ri'cnrrcnce  coefficients  7^  and 

<Tj  are  the  Schur  parameters  and  complementary  parameters,  respt  i  tively,  of  II.,  The  optimal 
polvnomial  is  then  given  by 

■  -I 

I', !>->{:)=  ^  f'u Ot_i(;). 

;.  =  i 
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where  the  vector  of  Fourier  coefficients  a  =  . Om-iF  ‘s  given  by  a  :=  (/*f  and 

f  [tni/(*l), ....  )]^-: 

Our  scheme  for  downdating  is  based  oii  the  observation  that  JI  and  0  can  be  computed 
from  II  and  U  by  applying  one  step  of  the  QR  algorithm  with  “ultimate”  shift  z.  together  with 
some  permutation  similarity  transformations,  to  determine  a  unitary  matri.x  IF  such  that 


’  i  0'^'  ]  [  ♦  croej  1  [  I  0’^  1  _ 

0  IF*  <roe,  Jl  J  [  0  IF  J  “ 


•  -T 

<Toe{  w 


(7oei  II 
w  0^ 


(2.10) 


U  0 


=  (/IF 


=  IF’a 


(2.11) 


Thus,  the  downdated  Fourier  cocfliciciits  are  obtained  by  accumulating  each  elementary 
unitary  transformation  against  a.  .Vn  elficieni  implementation  is  obtained  by  using  the  QR 
algorithm  for  unitary  llcssonberg  matrices  described  in  [■)]. 

We  now  assume  that  i  :*  for  some  .',!</<  m,  anil  let  tZ» :»  te/.  One  step  of  the  QR 
algorithm  with  shift  z  applied  to  the  matri.M  It  tloiermines  a  unitary  upper  Ilcsscnbcrg  matri.K 
such  that 

■  //  o' 

7/':= '/‘//r  = 


because  z  is  an  eigenvalue  of  II .  It  is  easy  to  si'o.  howovei',  that  the  QR  algorithm  applied  to  II 
will  not  yield  the  required  7/,  because  the  veil  or  fr^ei  will  not  be  transformed  ns  required  by 
(2.10).  On  the  other  hand,  an  RQ  algorithm  lor  77.  in  which  the  transforming  matri.x  V  would 
be  a  product  of  Givens  matrices  in  the  reverse  order  of  (2.12)  below,  would  yield  the  desired 
similarity  transformation. 

Instead  of  modifying  the  unitary  llcssonberg  QR  algorithm  of  [5]  to  perform  an  RQ  iteration, 
we  apply  the  QR  algorithm  to  the  unitary  Ib'ssenberg  matrix  //^  :=  .77/^7,  where  J  =s 
[eni,®m-i, — ®il  the  reversal  nialru  ol  orih'r  in.  |i  is  easily  seen  that  if  77  is  given  by  (2.5), 


77^  =  6'i(7i  )G’j(72)  •  ■ 

where  :=  ‘im-j'm  for  1  <  7  <  ur  und  :=  •,  The  application  of  iho  QR  algorithm  on 
77^  is  equivalent  with  that  of  the  RQ  nigoritiuii  on  11,  One  iteration  of  the  QR  algorithm  with 


shift  z  applied  to  11^  generates  a  unitary  upper  llesser.berg  matrix 


(2.12) 


such  that 


V'Il‘"V=  °  . 


(2.13) 


Moreover,  ^rn  be  taken  to  he  an  arbitrary  uniinodular  number  because  deflation  has  taken 
place  (o).  Let  tj  :=  (I  -  |«?j|*)’^'  denote  the  complementary  parameter  to  6j  for  1  <  j  <  m. 

Observe  that  only  the  last  two  components  of  l^Voem  are  nonzero,  and  that  |^m|  ^  1  can 
be  I'hosen  so  that  those  two  components  arc  given  by 


1  0 
I)  -^a 


^m-l  ^ni-l 
^111  — I  •'m— I 


0  _  r„_i(To 

m-  iko 


I'inally,  wo  transform  the  riglii-h.ami  ^idc  of  l>y  similarity  using  qY-  ^  ,  whore  J  is 

the  reversal  matrix  of  older  in  -  I,  imd  iranspo>e  the  result.  In  this  way  we  transform  Jl  to 
ir*//n'=  ir  ^  ,  whore  ir  =  ./r  i  °  ,  .Moreover,  IF VoCi  =  ^ 

0'  J  J  L  ®  ^  J 

by  the  uniqueness  of  the  reduction,  Jl”  *  //,  (To^'m-i  *  ^Ot  and  l^m-iko  *  te.  The  vector  of 
Fourier  coefficients  a  is  then  determined  from  (2.11)  by  applying  W*  to  a  incrementally. 

Observe  that  our  downdating  procedure  requires  knowledge  of  the  node  to  be  deleted  but  not 
of  the  corresponding  weight.  W'o  can  therefore  compare  the  computed  value  of  xb  to  the  actual 
value  wt  iit  order  to  assess  the  accuracy  of  the  computations.  If  w  is  not  close  to  wi,  then  the 
downdated  polynomials  6}  and  p,n-i  might  not  be  accurate.  .‘Vnother  accuracy  check  is  provided 
by  I  he  computed  v.due  of  /i,„  in  the  algorithm,  l  ids  (piantity  is  the  m'*‘  diagonal  element  of 
the  upper  triangular  matrix  in  the  QR  factori/atioii  of  II -  i/,  which  is  mathematically  zero 
when  i  is  an  eigenvalue  of  11.  If  the  compuiod  value  of  is  not  “tiny”,  then  the  computed 
downdated  polynomial  might  be  inaccurate.  We  therefore  consider  Pm  nnd  w  as  being  part  of 
the  output  of  the  algorithm. 


0 


Algorithm  2.1  (Downdating  by  removing  the  node-weight  pair  {i.w}  ) 
Input;  m.  {7j}7=ri> 

Output:  {<Tj}7=o’- 

h'jljLl*  •- 7m(7m-j]^l';  ! 

•“  [^»n-j+l]j=p' 

^0  ^0  •—  i» 

for  j  :=  -  I  do 

Pj  :=  {(f]  + 

if;  >  2  then  (7j_|  := 

’■/  (^j/Pj! 

fij  :=  ( +  7/j-i  )/Pj! 
nj  :=  -SjOj  + 

a,+i  :=  r,a,  + 

Aj  :=  l^j-i  +  I )//),:' 

jj  :=  Sjj  -  r;7j+i-7; 

Pm  +  7m^m-ll' 


^ni-l  !*  ^rn-lPm>  “  |^m-ll^0>'  ^0  * 

7m-;  :*  7m-i/|7m-i|;  <’m-i  :=  0;  (pnnunctor  coiroctioii) 


f,-,  1"'-;  im-l  . 


□ 


The  algorithm  overwrites  the  Fourier  coelficioius  {nj}’’L,  with  iuternirdiato  (luaiititios.  Al¬ 
gorithm  2.1  requires  0(»i)  arithmetic  operations  (  +  .  *,/)  and  tlie  evaliialion  of  ?n  -  1  square 

roots. 

We  liave  already  noted  that  if  the  nodes  arc  real,  tlien  the  analogue  of  the  niatri.\  II 
is  a  real  symmetric  tridiagonal  matrix  T  with  positive  subdiagonal  elements.  Similarly  with 
Algorithm  2.1,  downdatiiig  of  the  orthouormal  polynomials  associated  T,  as  well  as  of 
Fourier  cocfRcients,  can  be  carried  out  by  algorithms  based  on  the  QR  algorithm  for  symmetric 
iridiagonal  matrices.  This  observation  may  be  new., 

in  certain  data-filting  applications  it  may  be  desirable  to  update  the  polynomial  ;>m-i» 
given  by  ( l.;))-(1.10)  with  n  =  m.  by  replacing  certain  node-weight  pairs.  This  can  be  carried 
out  by  successively  removing  an  node-weight  pair  by  .Vlgorithm  2.1,  and  then  adding  a  new 
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iiode-woight  pair  using  one  step  of  Algorithm  •1.1  in  [j!].  This  combination  of  algoritliins  yields 
a  sliding  window  scheme.. 


3  Downdating  and  the  FFT  algorithm 

When  the  nodes  Zk  are  equidistant  or.  the  unit  circle  and  all  weights  trj  are  unity,  then 
interpolation  and  least-squares  approximation  of  a  given  function  by  algebraic  or  trigonometric 
polynomials  can  be  carried  out  rapidly  by  the  FFT  algorithm.  This  section  describes  in  some 
detail  how  our  downdating  scheme  may  be  combined  with  the  FFT  algorithm  to  achieve  rapid 
schemes  for  interpolation  when  the  nodes  arc  essentially  equidistant.  More  precisely,  we  will 
consider  the  case  when  the  set  of  nodes  {-fcli'Li  in  the  inner  product  (l.l)  is  a  subset  of  the 
set  of  e(|uidistaiit  nodes  / N ,  where  ><  :=  N  -  m  is  a  small  positive  integer..  In 

our  operation  count  we  will  a.'.sumc  that  n  is  iiidepeiideut  of  in.  The  weights  in  ( l.l)  are  all 
assumed  to  bo  unity. 

Let  /  dcuole  .1  loniph.'.x-N.ilueil  I'uik'Hoii,  v,  lio-e  value!-  /(;;,•),  1  <-  /•'  <  ni.  are  explicitly 
known.  We  remark  iliat  a  i('picsentaiion  of  ilie  interpolating  polynomial  //,„_)  €  Hm-i  in 
Newton  form  can  he  computed  in  O(in^)  arithmetic  opc'-ations.  The  Vandermonde  solver  by 
Djbrck  and  Poreyra  [2]  can  bo  used  to  determine  a  representation  of  p,n-i  in  terms  of  the 
monomial  basis,  and  requires  also  0(m^}  arithmetic  operations.  Our  scheme  only  requires 
0(m  log  m)  arithmetic  operations  and  yields  a  representation  of  pm-i  in  the  basis  of  Szego 
polynomials. 

Let  denote  the  complement  of  the'  set  in  {exp(2"(j/.V)}j\.~'.  and  let 

/(;[)  :=  0.  1  <  /r  <  K.  Thus.  /  is  defined  at  the  .V  lootsof  unity  oxp(27rd  J-  I  )/.V).  1  <  j  <  .V, 
and  we  can  compute  the  Fourier  coofficienls  of  the  polynomial  €  H.v-i  that  interpolates 
/  at  these  nodes  by  the  FFT  algorithm  in  0(.N’'  log  jV)  =  0(m  log  m)  arithmetic  oirerations. 
Using  the  Schur  parameters  given  in  Example  l.l.  we  apply  .-Mgorithm  2.)  k  limes  to  eliminate 
the  nodes  rj.,  \  <  k  <  k,  from  the  inner  product.  This  requires  0(m)  arithmodc  operations 
and  yields  the  Fourier  coellicients  of  the  desired  interpolating  polynomial  /),„-i.  The  Fourier 
coefficients  of  the  least  squares  approximants  p„-i  €  n„_i  with  n  <  m  are.  of  course,  a  subset 
of  the  Fourier  coellicients  of 

\  scheme  closely  related  to  the  one  outlined  .above  can  be  used  to  compute  trigonometric 
polynomials  rapidly.  Let  zi;  and  bo  the  noth's  introduced  above,  and  define  Oi,  :=  argfr/t), 
I  <  /.■  <  III.  and  :=  arg(c[).  1  <  A'  <  k.  .M.'O  .is.iiime  that  m  is  odd.  i.o..  in  =  2r  +  1.  Let 
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f(9)  be  a  real-valued  liiuctiou  defined  at  llic  notles  Ot,.  and  lot  /(^]t)  :=  0.  1  <  A;  <  k.  We  wish 
to  compute  a  trigonometric  polynomial  i(f))  €  ^pan{l.<:os  tf....,cos  jY^.sin  ^....,sin  r(9}  that 
minimizes  the  sum 

m 

This  can  be  accomplished  as  follows.  First  solve  the  minimization  problem 
min  f  1)^) 

by  the  FFT  algorithm  in  0(in  log  m)  arithnmtic  operations,  and  denote  the  optimal  coefficients 
by  ttj.  Remove  the  nodes  1  <  A:  <  k,  by  applying  the  downdating  scheme  k  times  to  the 
polynomial  pjr+K(-)  :=  where  A'  -  I  =  2r  +  k.  This  re(|iiires  0(m)  arithmetic 

operations  and  yields  the  polynomial  €  fljr-  I'he  desired  trigonometric  polynomial  is  then 
given  by  t(6)  :=  c“’  pi,.(c). ::  =  ('.\p(i^). 

4  Computed  examples 

We  now  present  the  results  of  some  numerical  e.\periments  with  our  downdating  procedure. 
These  e.xperiments  were  performed  in  FOllTR.V.V  on  a  SparcStation  SLC  at  Northern  Illinois 
University,  on  which  there  are  appro.ximately  7  and  10  significant  decimal  digits  in  single- 
precision  and  double-precision  calculations,  re.s  pec  lively. 

The  first  set  of  e.xperimcnls  compares  the  accuracy  of  the  downdating  procedure  with  that  of 
the  updating  procedure  lUQll  described  in  [S).  Wc  input  N  unimodular  nodes  positive 

weights  and  comple.x  function  values  For  any  positive  integer  m  <  N,  let 

am  denote  the  vector  of  Fourier  coefficients  of  the  .solution  Pm-i  of  (l-b)  with  n  =  m,  and  let 
gm  denote  the  vector  of  Schur  parameters  determined  by  the  inner  product  (1.1)  and  recurrence 
relations  {1.2)-(1.3). 

We  first  compute  vectors  a,v  and  gi\’  using  an  implementation  of  the  lUQR  algorithm  in 
single- precision  arithmeiic  [.S].  We  then  repeatedly  apply  a  single-precision  implementation  of 
Algorithm  2.1  to  compute  a,a  and  gm  for  decreasing  values  of  m.  The  A'*'  application  of  Al¬ 
gorithm  2.1  removes  llie  node  ;v-a+i  from  tin'  inner  |)in<lurt  to  compuio  the  solution  of  (l.(i) 
with  n  :=  m  :=  .V  -  A.  .After  each  dowmlating  step,  we  calculate  the  relative  error  in  dm,  i.e., 
||am  -  dmlli/llamlU,  and  the  error  in  gm,  i.e.,  |lg„.  -  gmlb,  where  ||x||2  denotes  the  Euclidean 
norm  of  x  €  C"*.  \\’o  also  solve  each  problem  of  order  m  using  the  lUQR  algorithm  in  single¬ 
precision  arithmetic  and  compute  the  resulting  errors.  The  results  of  the  lUQR  algorithm  in 


(Ion bio- precision  arithmetic  arc  used  as  exact  answers  in  the  error  calculations.  The  following 
tables  display  the  resulting  errors  for  the  dowiidaiing  procedure  (DD)  and  the  updating  proce¬ 
dure  ( UD).  In  all  of  the  cxpcriineuts,  each  function  value  f{zj)  has  its  real  part  and  imaginary 
part  randomly  generated  according  to  a  uniform  distribution  in  [-5,5]. 

We  first  choose  the  nodes  to  be  the  N  :=  300  roots  of  unity  Zj  :=  e.xp(2a’i(j  -  1)/N), 
1  $  ;  <  Jind  uniform  weights  wj  :=  1.  Table  1  shows  the  errors  for  problems  of  order 
in  ;=  iV  -  k  for  various  values  of  k.  Table  1  also  shows  the  results  with  the  same  choice  of  nodes 
and  weights,  e.xcept  that  the  nodes  are  permuted  in  a  random  way.  This  permutation  changes 
the  nodes  that  arc  deleted  as  well  as  the  order  in  which  the  nodes  arc  added  in  the  updating 
procedure.  It  should  Ix'  noted  that  the  eiror.s  in  the  downdating  procedure  can  be  expected  to 
increase  a.s  A'  increases.  Table  2  shows  that  similar  results  arc  obtained  with  the  same  set  of 
nodes  and  randomly  gcneratc<l  weights  in  (0,10). 

Table  3  shows  the  results  obtained  with  uniform  weights  and  .V  :=  300  nodes  zj  :=s 
exp(?r/(7  ~  l)/A')-  1  <  j  <  A  -  ill  [0,7r),  both  in  their  original  order  and  in  a  random  order. 
Here  again,  the  dowiulating  procedure  perfoini-  well. 

Table  4  shows  the  results  when  the  initial  .100  nodes  are  randomly  selected  points  on  the 
unit  circle.  In  this  example,  the  error  in  the  downdating  procedure  displays  a  sudden  increase 
at  A  s  10.  In  this  step  the  error  in  the  compiiKul  downdated  weight,  |u’  -  (h|,  was  greater  than 
10”'.  Observe  that  the  error  incurred  at  this  downdating  step  propagated  to  the  subsequent 
(lowiidating  steps,  but  the  errors  in  uiid  g,„  seem  to  grow  gradually  ns  A  increases,  i.c.,  as 
in  N  -  A  docrea.se.s.  Other  experiments  with  random  nodes  produced  similar  results.  It 
should  ho  noted  that  in  our  cxporimciits.  n  Inige  error  in  the  computed  weight  did  not  always 
coincide  with  a  large  Jump  in  ilic  errors  in  a,„  and  g,„.  It  is  clear  that  more  work  needs  to  be 
done  in  order  to  understand  the  numerical  aspects  of  the  updating  and  downdating  problems. 

Our  final  experiment  tests  tiie  accuracy  and  speed  of  the  procedure  described  in  Section 
3  for  downdating  the  l^FT.  The  .'V-point  FFT  is  used  to  obtain  the  Fourier  coefficients  of  the 
interpolating  polynomial  P;v-i.  -‘V  raiulondy  selected  set  of  nodes  is  then  removed  from  the 
inner  product  using  .Algorithm  2.1.  'I'liis  ( xpeiiiiiont  was  run  with  a  radi.x-2  FFT  subroutine  in 
single  precision  iirilhnioiic  with  .V  :=  1021  and  .V  ;=  2018.  Table  5  shows  the  computed  error 
after  A  downdating  steps  for  various  values  of  A.,  As  above,  wo  use  the  results  of  the  lUQR 
algorithm  in  double  precision  aiitlimelic  as  ox.ict  answers  for  error  checking.  We  also  display 
the  amount  of  CPU  time  required  by  the  FFT  with  A  downdating  steps  and  the  lime  required 
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Table  1:  300  nodes  with  equispaced  arguments  in  [0,2ff);  uniform  weights. 


nodes  in  order  of  increasing  argument. 

11 - 

nodes  in  random  order. 

relative  error  in  a,„ 

relative  error  in  u,„ 

k 

l)D 

ID 

WmSm 

DD 

UD 

0 

U.i22t:-03 

0.122E-03 

0.12()E-0.3 

U,l2()E-03 

0.2S4E-04 

0.284E-04 

0.310E-04 

0.340E-04 

10 

0.141E-03 

0.121E-03 

().29.5E.03 

u.iior,-03 

0.333E-04 

0.290E-04 

0.545E-04 

0.340E-04 

20 

O.MlE-03 

0,115E.03 

0.405E-03 

0.118E.03 

0.473E-04 

0.314E-04 

0.211E-03 

0.361E-04 

30 

0.11.5E-03 

O.lOOE-03 

U.322E-03 

0.1o4E-03 

0.550E-04 

0.272E-04 

0.272E-03 

0.331E-04 

40 

0.112E-03 

0,9(57E-04 

0.340E-03 

0.1.55E-03 

0.943E-04 

0.257E.04 

0.359E-03 

0.374E-04 

50 

0.112E-03 

0.940E-04 

0,320E-03 

0.164E-03 

0.964E.04 

0.259E-04 

0.325E-03 

0.437E.04 

70 

0.124E-03 

0.102E-03 

0.578E-03 

0.171E.03 

0.174E-03 

0.258E-04 

0.424E-03 

0.422E*04 

90 

0.132E-03 

0.10.')E-03 

0.725E-03 

0.173E-03 

0.254E-03 

0.235E-04 

0.677E-03 

0.422E-04 

110 

0.154E-03 

0.107E-03 

0.365E-03 

0,104E-03 

0.375E-03 

0.238E-04 

0.753E-03 

0.383E-04 

130 

0.151E-03 

0.U7E-03 

0.289E-03 

O.lGOE-03 

0.465E-03 

0.226E-04 

0.109E-02 

0.3.59E.04 

1.50 

O.lOCE-03 

0.127E-03 

0..‘121E-03 

0.1.50E-03 

0,600E-03 

0.186E-04 

0,119E.02 

0.4G3E.04 

by  the  single-precision  IU(Jll  algorithm  on  the  i>rol)lcin  of  order  m  =  .V  -  k,.  It  is  interesting  to 
note  tliat  the  dowiidaiiiig  procodnre  produces  siihsianlially  more  accurate  answers  faster  tliau 
tlic  lUQR  algorithm  evem  for  moderately  si/cd  values  o{  k.: 

5  Conclusion 

New  algorithms  are  presented  for  discrete  least  squares  approximation  by  algebraic  poly¬ 
nomials  on  the  unit  circle  and  by  trigonometric  polynomials.  The  algoritluns  downdaie  the 
•  approximant  as  well  as  recurrence  cocfTicicnts  for  Szego  polynomials  when  a  node-abscissa  pair 
is  removed  from  tlie  inner  product.  The  schemes  are  based  on  the  QR  algorithm  for  unitary 
llesscubcrg  matrices.  particularly  attractive  application  is  to  combine  our  algorithm  for 
trigonometric  polynomials  with  the  fast  rourier  transform  algorithm.  This  yields  a  fast  scheme 
Cor  computing  trigonometric  least  squares  approximants  when  the  nodes  of  the  inner  product 
form  a  subset  of  c((uidistant  nodes  on  (0,  2”). 
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Tabic  2:  300  nodes  with  equispaced  arguments  in  [0,2;r);  random  weights  in  (0,10) 


nodes  in  order  of  increasing  argument. 

nodes  in  random  order. 

relative  error  in  u,n 

error  in  g,„  I 

relative  error  in  Am 

k 

DD 

UD 

DD 

UD 

DD  , 

UD 

0 

0.207E-03 

U.207E-03 

0.2'ME.03 

0.211  E-03 

0.511E-04 

0.511E-U4 

0.572E-04 

0.372E-0-1 

10 

0.200  E-03 

0.1S3E-03 

0.33.)E-03 

0.181  E-03 

0.363E-04 

0.537E.04 

0.844E-04 

0.744E-OI 

20 

0.199E-03 

0.171  E-03 

O.‘!17E-03 

0.192E-03 

0.C29E-04 

0.512E-04 

0.140E-03 

0.722E.0-I 

30 

0.187E-03 

0.175E-03 

0.302  E-03 

0.198E-03 

0.78oE.01 

0.518E-04 

0.252E-03 

0.750E-04 

•10 

0.188E-03 

0.1 79  E-03 

0.3'l0E-03 

0.203E-03 

0.873E-04 

0.509E-04 

0.303C-03 

0.739E.04 

.*)0 

0.181  E-03 

0.107E-03 

().33'lE-03 

0.209E-03 

0.308E.03 

0.487E-01 

0.4 14  E-03 

0.714E-0-I 

70 

0.208  E-03 

0.1()9E.03 

0,(i0-lE-03 

0.217E.03 

0.406E-03 

0.479E.04 

0.C21E-03 

0.G38E.0-1 

90 

0.192E-03 

0  l  l.')E-U3 

0.70.3  E-03 

0.222E-03 

0.572E.03 

0.442E-04 

0.178E-02 

0.688E-0t 

110 

0.1'J8E-03 

U.15oE-03 

0.389  E-03 

0.220E-03 

0.874E.03 

0.333E-04 

0.184E-02 

0.533E-0'I 

130 

0.2U1E-03 

0.155E-03 

0.341E-03 

0.213E-03 

0.694E-03 

0.286E-04 

0.139E-02 

0.401E.04 

130 

0.201E-03 

O.lOOE-03 

0,379E-03 

0.195E.03 

0,122E-02 

0.256E-04 

0.153E-02 

0,285E-04 

Table  3:;  300  nodes  with  equispaced  arguments  in  [0,  tr);  uniform  wciglils. 


nodes  in  order  of  increasing  argument. 

nodes  in  random  order. 

lelativc  error  in  n,„ 

relative  error  in  a,n 

k 

DD 

UD 

iiyi 

DD 

UD 

■Elaii 

0 

0.398E-03 

0.398E-03 

0.712E-03 

0.742E-03 

0.125E-03 

0.125E-03 

O.llOE-03 

0.116E-03 

10 

0.455E-03 

0.398E-03 

0.188E-02 

0.7.13E-03 

0.120E.02 

O.lOlE-03 

O.lOCE-02 

0.124E-0.3 

20 

0.437E-03 

0.389E-03 

0.187E-02 

0.720E-03 

0.M4E-02 

0.112E-03 

0.146E-02 

0.122E-0.3 

30 

0.452E-03 

0.390E-03 

0.228E-02 

0,7 14  E-03 

0.170E-02 

O.llOE-03 

0.163E-02 

0.131E.03 

40 

0.479E-03 

0.391  E-03 

0.265E-02 

0.G94E-03 

0.195E-02 

0.947E-04 

0.107E-02 

0.U3E-0.3 

.50 

0.508E-03 

0.391  E-03 

0.271  E-02 

0,080E-03 

0.215E-02 

0.886E-01 

0.105E-02 

0.109E-0.3 

70 

0.57CE-03 

0  131  E-03 

0.27SE-02 

0.04.3E-03 

0.251E-02 

0.796E-04 

0.210E-02 

0.109  E-03 

90 

0.533E-03 

0.450E-03 

0,1 08  E-02 

0,00.3E-0.3 

0.307E.02 

0.849E-04 

0.302E-02 

0.106E.0.3 

no 

0.558E-03 

l).13.3E-03 

0  193  E-02 

0..3.)2E03 

0.4.38E-02 

0,795E-04 

0.307E-02 

0.102E-I)3 

1.30 

0.f)50E-03 

0  l23i:-()3 

0  •JIOE-O'J 

0  303 E-o;'. 

0..32SE-02 

0.792E-0 1 

0.339E-02 

U.840E-0I 

1.30 

O.llOE-02 

0  lO.lE-Od 

0  .3.30  E-02 

0  l()0E-03 

0.071  E-02 

0.548E-01 

0.528  E-02 

().895E-01 
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Tabic  •!::  300  nodes  with  raiidoia  armimciiis  in  [0,2Jr);  uniform  weights. 


k 

I'l'lative  error  iii  a,„ 
l)D  fD 

error 

DD 

ii‘  «m 

UD 

n 

i)  iiii:.03 

0.  lllE-t)3 

0  yj5E-03 

0.y25E-03 

0  .j0j!C.03 

U.  10Ui:-03 

0.9ME-03 

0.838E-03 

10 

0  JlJiE-Ol 

O.lGOE-03 

O.aiSE-Ol 

0.838E-03 

•JO 

0  '.njE-oi 

0  rjOE-03 

0.373E-01 

0.7S0E-03 

30 

0.'2'>OE-01 

0  1(5‘IE.03 

O.3’28E.0l 

0.735E-03 

•to 

0  lllE-01 

0.117I>03 

0.131E-01 

0.826E-03 

•jO 

O.IJJOE.OI 

O.r.ME-03 

0.2‘10E.01 

0,G54E-03 

70 

0  lOOE-Ol 

0.137E-03 

0.’2e<lE.01 

0.658E-03 

90 

0  JIOE-Ol 

0.1’J’JE.03 

0.'275E.0l 

0.483E-03 

110 

0.317E.01 

0.137E.03 

0.8-20E.01 

O.IOOE.02 

130 

O.50-2E-01 

0,r20E-03 

0.85JE-01 

0.508E-03 

150 

0.513E.01 

0.7SSE-0-1 

0.8G5E-01 

0.3G0E-03 

V 


Table  5:  Downdating  the  FFT 


.V  = 

1024 

CPU  seconds  I 

lU 

k 

DD 

UD 

IQI 

11 

U.-IOIE-Oo 

U.213E-01 

0.218E-05 

0.204E-01 

0.303E+01 

0.M8E+03 

■i 

0.502E-05 

0,215E-01 

0.623E-05 

0.195E-01 

0.814E-I-01 

0.142E-i-03 

E 1 

■1 

O.fJOSE-O') 

0.217E-01 

0.M2E-O4 

0.103E-01 

0.132E<l-02 

0.136E+03 

923 

101 

O.yOOE-05 

0.l93E-0i 

0.257E-04 

0.244E-01 

0.253E-f02 

0.122E+03 

913 

111 

0.9'l<tE-0r> 

0.191E-01 

0.205E-04 

0.242E-01 

0.276E+02 

0.119E+03 

893 

131 

O.lllE-Oi 

0.178E-01 

0.312E-04 

0.219E-01 

0.323E+02 

0.114E+03 

873 

151 

0,129E-0'l 

0.17-lE-Ol 

0.373E.04 

0.220E-01 

0.368Et02 

0.109E-f03 

853 

171 

O.MOE-O'l 

0.172E-01 

0.462E.04 

0.221E-01 

0.412E-f02 

0.104E+03 

83w 

191 

0,172E-01 

0.17-1E-01 

0.566E-04 

0.220E-01 

0.456E+02 

0.990E+02 

S13 

211 

0.223E-01 

0.171E.01 

0.721E-04 

0.241E-01 

0.407C-f02 

0,941E+02 

713 

311 

1).  1771- -01 

(),133E.01 

0.131E-03 

0.169E*01 

0.692E+02 

0.719E+02 

(ii3 

•111 

0.107E-():i 

O.lO-lE-Ol 

O.I5SE-03 

0.275E-01 

0.862E+02 

0.527E-f02 

513 

.511 

().2C7E-()3 

().709E-02 

0.527E.03 

0.174E-01 

O.lOlE+03 

0.365E+02 

113 

(ill 

lJ.28(iE-03 

0.117E-02 

0.I87E-03 

0.852E-02 

0.112E-I-03 

0.233E+02 

313 

711 

0.337E-03 

0.322E02 

0.554  E.03 

0.781E.02 

0.122E+03 

0.132E+02 

leliitive  error  in  iin, 

ilH 

m 

DD 

UD 

2037 

11 

U.598E-05 

0.105E+U0 

rnsmm 

ItlT:n3gSiril 

0.391E+01 

0.390E+03 

1997 

51 

0.7.I0E-05 

().ll()E+00 

0,151E.04 

0.235E+00 

0.172E+02 

0.374E+03 

1947 

101 

O.IOOE-Ol 

0.102E+00 

0.293E-04 

0.164E+00 

0.33CE+02 

0.355E+03 

1897 

151 

0.125E-Ut 

().y80E-01 

0.420E-04 

0.126E+00 

0.495E+02 

0.337E+03 

1847 

().159E-0I 

O.OISE-Ol 

0  O-PiE-O  l 

0.131E+00 

0.649E+02 

0.319E+03 

1797 

251 

().237|.-n.| 

0.881E-01 

().8I2E.04 

0.136E+00 

0.799E+02 

0.301E+03 

1747 

301 

0  2>!)lv()l 

0  805C-()1 

0  1  lOE-03 

0.157E+00 

0.9I6E+02 

0.285E+03 

1097 

351 

1),;158I>0I 

1)..S31E-01 

().ll5E-03 

0.137E+00 

0.109E+03 

0.268E+03 

1647 

•101 

U.51.|E-()1 

0  fiOOI'-Ol 

0  170E.03 

0.125E+00 

0.123E+03 

0.252E+03 

1.597 

151 

0  749E-()  1 

0  748li:-01 

0.242E-03 

0.133E+00 

0.13CE+03 

0.237E+03 

1.547 

501 

0  779E-04 

l).729E-01 

0.243E-03 

0.132E+00 

0.149E+03 

0.222E+03 

1447 

001 

0  953E.()  1 

0.095C-01 

0  313E.03 

U.149E4-00 

0.174E+03 

0.194E+03 

1347 

701 

0.242E-03 

0 028E-01 

0.396 E-03 

0.114E+00 

0.197E+03 

0.1'j7E+03 

1147 

901 

0  .517E-03 

0.443E.01 

0.109E-02 

0.116E+00 

0.238E+03 

0.i:iE+03 

1047 

1001 

U.73lE-()3 

0.381E-01 

0.113E.02 

0.640E-01 

0.256E+03 

O.luOE+03 

947 

1101 

0.114E-()2 

0.285E.01 

0.183E.02 

0.523E-01 

0.273E+03 

0.815E+02 

\  >• 
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